Decomposition of modified Landau-Lifshitz-Gilbert equation 
and corresponding analytic solutions 
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The Suzuki- Trotter decomposition in general allows one to divide the equation of motion of a 
dynamical system into smaller parts whose integration are easier than the original equation. In 
this study, we first rewrite by employing feasible approximations the modified Landau-Lifshitz- 
Gilbert equation for localized spins in a suitable form for simulations using the Suzuki- Trotter 
decomposition. Next we decompose the equation into parts and demonstrate that the parts are 
classified into three groups, each of which can be solved exactly. Since the modified Landau-Lifshitz- 
Gilbert equation from which we start is in rather a general form, simulations of spin dynamics in 
various systems accompanying only small numerical errors are possible. 
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I. INTRODUCTION 

The recent development of fabrication techniques of 
nanostmctures and the extension of knowledge about 
magnetic properties of condensed matters are accelerat- 
ing the efforts for achieving more efficient electric control 
of magnetic states for applications in engineering. The 
Landau-Lifshitz-Gilbert (LLG) equation[l], Q has been 
widely used for interpretation of experiments and theo- 
retical analyses for complex magnetization dynamics. In 
particular, its phenomenologically modified versions 0-[H 
proposed for incorporating the effects of spin torque in- 
duced by an electric current |(| are important both for 
realization of spintronics devices and for micromagnetic 
simulations of such nonequilibrium phenomena. Theo- 
retical studies, including numerical simulations, often ad- 
dress the space in which magnetic moments reside as a 
continuum since the typical length scale of magnetization 
patterns is on the order of ten nanometers, much larger 
than lattice constants. 

On the other hand, there are many systems whose spin 
configurations are known to be crucially influenced by 
the microscopic arrangement of each spin such as non- 
collinear spin systems, which include helical magnetism, 
weak ferromagnetism, and frustrated spins. The realiza- 
tion of spin configuration in such systems often originate 
from the symmetry of crystal lattice and some of them 
are related to relativistic effects on electrons, which are 
in general effective only in the vicinity of each ion, de- 
manding microscopic considerations. 

For construction of a reliable model for a specific elec- 
tronic system, schemes for determination of its param- 
eters using first-principles electronic structure calcula- 
tions have been developed. They allows one to calculate 
the effective exchange integrals [7| and the Dzyaloshinskii- 
Moriya interaction |8| between localized spins. Schemes 
taking the electronic correlation effects into account for 
calculations of the Coulomb repulsion and the Coulomb 
exchange integrals @ between localized electronic orbitals 
have also been developed. Even for spin simulations of 
large length scales, those which are the most reliable 



should be performed using model parameters determined 
by first-principle calculations. To perform simulations 
of spin dynamics by incorporating the effects coming 
from the microscopic arrangement of spins and evaluate 
quantitatively the properties of the system, not only the 
material-specific parameters but also reliable algorithms 
for numerical integration are needed. One of the most re- 
liable integration scheme is the Suzuki- Trotter decompo- 
sition (STD) method^, whose numerical error is third- 
order with respect to the time step. If we know the exact 
solutions of the individual parts of a decomposed time- 
development operator, the STD method ensures that the 
source of numerical errors is only the STD method itself. 

This paper is organized as follows. In Sec. II, we 
rewrite the modified LLG equation into a physically ap- 
propriate form for localized spins, which is demonstrated 
to consist of terms of three kinds. In Sec. Ill, the an- 
alytic solution of the equation of motion for each of the 
parts is provided. 



II. MODIFIED LANDAU-LIFSHITZ-GILBERT 
EQUATION 

The LLG equation and its modified versions are for 
magnetization M, which is defined as the total magnetic 
moment density as a function of continuous space coor- 
dinate. In the present study, we use the effective spin 
angular momentum S = —vM/j, where v is an appro- 
priate volume unit, rather than M. v depends on the 
spatial scale for physical phenomena one would like to 
see and study. We assume that the effective spin be- 
haves as a single classical spin. Although the validity 
of the choice of v and the treatment of the single spin 
should be of course confirmed for individual simulations 
of specific systems, such problems are out of scope of the 
present study. 

In this subsection, we rewrite the phenomenologically 
modified LLG equation by replacing the spatial deriva- 
tives with finite differences of spin variables for two rea- 
sons. The first is that we are treating the system as 
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consisting of localized spins which are not necessarily ar- 
ranged on a regular mesh. The second is that the rewrit- 
ten parts of the equation of motion are of the form whose 
exact solutions are known. 



A. Original form 

We express the effective external field felt by the n-th 
spin in the system as H n = ^B n - J2 m =i n JnmS m - 7 is 
the gyromagnetic ratio and B n is the external magnetic 
field. Jnm is a 3 x 3 matrix representing the effective 
exchange interaction between the n-th and m-th spins, 
m in the summation runs over all the spins surrounding 
the n-th spin. By using a symmetric 3x3 matrix J7„ , we 
express the field due to the magnetocrystalline (single- 
ion) anisotropy as ~2J^S n: whose factor 2 comes from 
the quadratic form of the spin Hamiltonian. [ll], [l2[ We 
introduce constraints on the magnitudes of the spins in 
the present study as S n = const. Thus can be as- 
sumed to be traceless without loss of generality. 

Kim et al. Q recently derived a general form of the 
LLG equation for studying magnetization dynamics in 
a thin ferromagnet incorporating the Rashba effect [l3j|. 
which is a kind of spin-orbit interaction due to lack of in- 
version symmetry. Meanwhile, Zhang and ZhangQ had 
proposed an enhanced damping tensor, which gives rise 
to temporarily and spatially varying damping term, as 
a modification to the LLG equation. It was introduced 
to describe the effect of a conducting electron in a fer- 
romagnet which carries away the excess angular momen- 
tum of the precessing ferromagnetic moment. [l4l. [l5j The 
current-driven rotational motion of a domain wall is in- 
fluenced significantly by the damping enhancement . (l6| 

We adopt the following form of the modified LLG equa- 
tion in the most general form to date as a starting point 
of the present work: 



^ = S n X \ H n + if R + ^-S n X H R 



-2S n X J*S n + (u • V)5„ + ^S n x (u • V)S 



1 



Sn 
dS n 



, Sn X V 

o n at 



(1) 



where 



Vu = aSij + ^ yVS„ x d k S n ) t (S n x d k Sn)j (2) 
o n 



k = x,y, z) is the enhanced damping tensorQ with 
the dimcnsionlcss Gilbert damping constant a and the 
parameter n in a dimension of squared length. The elec- 
tric current density j e and its polarization rate P defines 
the velocity vector u = n&P/[M s {\ + P 2 )]j e , where (3 is 
the ratio of the nonadiabatic spin transfer torque to the 
adiabatic one0, [l7| and M s is the saturation magneti- 
zation, ub is the Bohr magneton. H R = —k^{e z x u) 



is the effective field introduced for the Rashba effect, Q 
whose strength is measured by /cr in a dimension of in- 
verse length. 

Since we have adopted the expression of the effective 
field as the summation of the contributions from the sur- 
rounding effective spins, our formulation allows for an ar- 
bitrary mesh of space for coarse-graining. Furthermore, 
the derivation of the equation of motion below is appli- 
cable to microscopic (not effective) spins distributed in a 
crystal. 



B. Rewritten form 

As stated above, the spatial derivatives in Eq. ([1]) 
should be replaced with finite differences appropriate for 
specific systems. One might simply calculate the spa- 
tial derivative of a spin direction by using its neighbor- 



a<s. 



n • 

nm 



ing ones as diS n ~ AfS n + 

The appropriate definitions of the vectors A™ and A 
which characterize the numerical derivative operator Aj 
depend on the locations of the spins in the system, in- 
cluding the dimensionality and the periodicity. The 
constant-magnitude condition of the spin requires that 
the spin and its spatial derivative be orthogonal since 
S n ■ diS n = diS n /2 = 0. The numerical derivative pro- 
posed above, however, does not ensure the orthogonality. 
We therefore subtract the component parallel to S„ from 
A t S n as Aj-Sn = A t S n - (S n ■ A l 5„)5„/^ 2 i . We can 
thus write 



where 



AlS n = A l S n -A"S n , 



AiS n = K m S r , 
X\\c _ a nm S n- S m 



(3) 

(4) 
(5) 



It is noted that AiS n does not contain the variable S n . 
We adopt the constant-magnitude numerical derivative 
Aj-Sn instead of d{S n in the present study. 

Let us then rewrite the equation of motion by employ- 
ing the constant-magnitude numerical derivative. There 

clearly exists a relation Aj-S„ = —S n x (S n x AiS n )/S n , 
which can be used for rewriting the spin transfer torque 
terms in Eq. ([T]) as 



(« • V)S„ + -|-5 n x (« • V)S n 

O n 

= PS n x H n ^~S n x (S n x H n )i 

On 



where we have defined 
1 



O n 



(6) 



(7) 
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The equation of motion of the form Eq. ([T} is not 
tractable since the time derivative appears on both sides. 
To make the equation computationally manageable, we 
operate S n x T> on both sides of Eq. ((TJ). We calculate 
the right hand side below by dividing it into four parts. 

The first and second parts are calculated as follows. 
We neglect the terms which involves r]H R and r]J° since 
they originate from the spin-orbit interactions and their 
energy scale is much smaller than those for nonrelativis- 
tic interactions in general and the effects of these terms 
on the qualitative behavior of enhanced damping are ex- 
pected to be insignificant. Thus we can rewrite the first 
part as 



S„ x V 



S n x lH n +H R + -^-S n x H R 



= aS„ x [S n x (H„ + H R )} - af3S n S n x H R 

+^J2(^Sn ■ H n )S n X (S„ X A t S n ) (8) 

n i 

and the second part as — 2S n x T>(S n x J^Sn) = 
— 2aS n x (S n x JnS n ). Furthermore, we neglect terms 
which involves f3rj for the third part since /3 is much 
smaller than unity ljj , and it is rewritten as 



S n x V(u • V)S„ + -J-S„ x V[S n x (u • V)S„] 

On 

= aS n S n x H n + af3S n x (S n x H n ) 

+ lkY. Sn - x Hn)S n X (S n X AiS n ). (9) 

n i 

The fourth part to the first order of rj is rewritten as 



a S; 



ns n (Ais n )(Ais n y 




(10) 



Combining these four parts, we obtain 

S n X 



dt 



= aS n X [S n X (H n + H R )] 

-aPSnSn x H R 

+ ^J2(AiS n ■ H n )S n X (S n X AiS n ) 
?>n i 

~2aS n X (S n X J^Sn) 

X H n + aftSn X (S n X H n ) 
""^J ^ ' S n • (AiS„ X H n )S n x (S„ X AiS n ) 



-a 2 S n 



dS n 

dt 



(11) 



We have omitted the term involving 
(ar]/S 2 ) J+ \A^-S n \ 2 coupled to the time derivative 
since this term is at most comparable to a 2 [3, EH , which 
is much smaller than unity, and it would not affect 
significantly the entire time derivative. Decomposing 
A^S n as Eq. ([3|) and H n into components parallel 

and perpendicular to AiS n as H n = H ni + H^ t , 
we can write, using Eq. JS]), (A^S n ■ H n ) = 

(AiS n ■ H n ) - A]„(S n • H^) - Aj n (Sn ■ H^). 
By defining the two kinds of functionals as 
A [S;C,a,b] = aS x C + bS x (S x C) and 
Ai[S;C,6] = (b ■ S)S x (S x C)/S , and substi- 
tuting Eq. pT|) back into the original equation, Eq. ((T|), 
we finally reach the following form: 



dS_ 

~dt 



^(eff) + ^(MCA) + ^(SPC) + A (SOI) 

+ A (ED0) +A (ED| ! ) +A (ED±) j (12) 



where 



A (off) = A [S;H,K,-Ka/S], 



(13) 



A (MCA) = -2kS x J S S + 2^-S x (S x J S S), (14) 

A (SPC) = Aq [ S . jj K (_ a + ^ _ K ( a/3 + 1 ^ s j j ( 15 ) 

A (SOI) = Aq [ S . Hr k(1 + + (lg) 

a (edo) = A Q [S; 0, -Kr)(AiS ■ H)/S% (17) 



a™ = J2MS;\s, kv a!h!/s% (is) 



(ED_L 



(19) 

and k = 1/(1 + a 2 ). We have omitted the subscript 
n. The equation of motion obtained here consists of 13 
parts, which are only of three kinds, Aq, A\, and the 
magnetocrystalline anisotropy terms. This is the first 
main result of the present work. Aq describes an ordinary 
precessing and damped motion of a spin. A\ describes 
an abnormal damped motion, which cannot be expressed 
as a special case of Aq. Since rj is no longer coupled 
to the time derivative, the damping enhancement can be 
treated on the same footing as the ordinary precessing 
and damped motion. 



III. EXACT SOLUTIONS 

We have obtained the equation of motion for a generic 
case, Eq. (fT2)l , however, we could not hope that its ex- 
act solution for arbitrary parameters are obtained. The 
equation should thus be solved numerically. Ma and 
Dudarev[H| adopted the STD method [13] for the nu- 
merical integration of the equation of motion of spins. 
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A first-order differential equation of a dynamical vector 
variable x, 

— = (L 1+ L 2 )x, (20) 

where L\ and L 2 are operators acting on x, has a formal 
solution x(t + At) = exp[(Li + L 2 )At]x(t). The STD 
method allows one to decompose the time evolution op- 
erator as 

e (L 1+ L 2 )At = ^ At/2 e L 2 At ^ At/2 + ( At 3) j ( 2 1) 

which means that if we divide a complex equation of mo- 
tion into sufficiently small parts whose exact solutions 
are known, the numerical errors of a simulation of the 
dynamics are restricted within those due to the STD 
method itself, as small as third-order with respect to the 
time step. The advantages of the STD method over other 
integration methods for spin dynamics are explained in 
detail in literature. [HI In what follows, the exact solu- 
tions of the three kinds of equations of motion derived 
above are provided. While the exact solution for A has 
been reported, those for A\ and the magnetocrystallinc 
anisotropy term with general parameters are, to my best 
knowledge, first to be provided as the second main result 
of the present work. 



as a function of t, we have to examine two cases accord- 
ing to whether b is parallel or perpendicular to C. Even 
when b is neither parallel nor perpendicular to C, decom- 
position of b into the two components allows us to use 
the STD method, which is the reason for the separation 
of A™ and A< ED± ). 

1. For b parallel to C 

When b is parallel to C, the equation 9 must satisfy is 
c\f) 

—— = Sb z C cos 9 sin 9. (25) 
at v ' 

By integrating both sides, wc obtain 9 = 
arctan[tan#o exp(Sb z Ct)], where 6q = 9(t = 0). If 
b z C > 0, 9 = 7r/2 is a stable stationary point and 
9 = and it are unstable stationary points. If b z C < 0, 
6 = 7r/2 is an unstable stationary point and 9 = and tt 
are stable stationary points. (See Fig. [T]). 

The solution for C in an arbitrary direction is obtained 
by rotating correctly the solution given above. Its expres- 
sion is 

a _S X (e-«-l)C + CS 
CWx 2 e~ 2Ct + 1 - X 2 
where C = S(b ■ C) and X=(S - C)/{SC). 



A. Ordinary precession and damping term 



2. For b perpendicular to C 



Since the exact solution of the equation of motion gov- 
erned only by an Aq term, 



dS_ 



aSxC + bSx{Sx C), 



(22) 



has been provided in an earlier work[l8f . only its expres- 
sion is written here. With the parameters Sq = S(t = 
0), £ = aC, C = bCS, X=(S - C)/(SC), it is given by 



C[l + e 2 C 4 + x(i-e 2 C*)] 
■[20e ct cos^So + 2e c * sin£i(S x C) 
+S{1 - e 2<t + x(l + e 2Ct - 2e c *cosf*)}C]. (23) 

B. Abnormal damping term 

Let us consider the equation of motion governed only 
by an A\ term, 



(24) 



We first assume C to be along the z-axis: C = Ce z . 
Expressing the equation of motion with the spherical co- 
ordinates 9 and qb for S, we immediately obtain <fi = 
const. = 0o- F° r obtaining the analytic expression of 9 



When b is perpendicular to C, the equation 9 must 
satisfy is 



4- = 5fe ± Csin 2 i 
at 



(27) 



where b± = b x cos(f)o + b y shi(f>Q. By integrating both 
sides, we obtain 9 = arctan[(cot 9q — Sb^Ct)" 1 }. If 
b±C > 0, 9 = is an unstable stationary point and 
9 = 7r is a stable stationary point. If b±C < 0, 9 = is a 
stable stationary point and 9 = 7r is an unstable station- 
ary point. It is noted that 9 = tt/2 is not a stationary 
point. (See Fig. p. 

The solution for C in an arbitrary direction is obtained 
by rotating correctly the solution given above. Its expres- 
sion is 

-sttc + cs 

cv(x-a) 2 + i-x 2 

where ( = C(b- S ) and X =(S - C)/{SC). 

C. Magnetocrystalline anisotropy term 

Let us consider the equation of motion governed only 
by a magnetocrystallinc anisotropy term, Eq. (|14[) . 

? = -2kS x J s S + 2-^S x {S x J S S), (29) 
at b z 
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(a) .4o 
&>0 kC 



b<0 kC 



(b)Mb\\c) 
b-c>Q kC 




Eq. ([29)1 becomes via operation of R on it from left as 



(c) Ax (61c) 



FIG. 1: (Color online) Directions of ordinary and abnormal 
damped motion of S for equations of motion governed by (a) 
Ao and (b), (c) A\. Green (brighter) arrows represent spins 
at unstationary points, while blue (darker) arrows represent 
stable or unstable stationary points. In (a), precession is as- 
sumed to be absent (a = 0). In (c), b, C and S are assumed 
to lie on the same plane. 



— ^ = — 2t7(A 3 - X 2 )S y S z 
"2-J^[AiS' 2 — {\\Sl + \-2Sy + \3S^)]S X , 

1 Q 

= -2r?(A! - \ 3 )S Z S X 

~^^2 i^ 2 ^ 2 _ (^l^x + ^Sy + X^S^ )]S y , 



dS 1 

— -i = -27y(A 2 - Xi)S x S y 



~^g2 i^S 2 — (Ai-S 1 ^ + X2S 2 + \3Sz)]S z 



(31) 



(32) 



(33) 



where S = RS and r\ = k det R. Whether this system of 
equations can be solved exactly depends on the eigenval- 
ues. If the three are equal, S does not change. If two of 
them are equal and the other is different, the equations 
of motion, Eqs. (|3"Tj) - (|33|) . can be solved exactly. If the 
three arc different from each other, exact solutions can- 
not be obtained and we resort to dividing the equation of 
motion into two parts, the precession and the damping 
terms, whose exact solutions can be obtained separately 
Thanks to the STD method, this division only increase 
the number of time-development operators and not acts 
as a source of numerical errors. Thus the consideration 
for the following three cases suffices. 



1. For Ai = A2 



where £ = naS. This equation describes the dynamics 
of a spin feeling a torque coming from its own direction. 
Since J s is a symmetric matrix, it can be diagonalized 
by a real orthogonal matrix R. Let Vi(i = 1,2,3) be 
the eigenvector and Ai be the corresponding eigenvalue. 
Then we can write *i? = (vi, i?2, U3) and J s = t RDR, 
where D = diag[Ai, A2, A3]. *i? is the transposition of 
R. If the eigenvalues Ai, A2, and A3 are different from 
each other, we require that the eigenvalues be arranged 
according to the initial spin direction so that 



Ai > A 2 > A 3 (S ■ J S S Q - X 2 S 2 < 0) 
Ai < A 2 < A 3 (So ■ J S S Q - X2S 2 > 0) 



(30) 



for later convenience for the £ = case discussed below. 



Here we consider a case in which two of the eigenval- 
ues are equal. We can assume that Ai = A 2 = A and 
det R = 1 without loss of generality. We use the spheri- 
cal coordinates 9 and cj> of S. The equation 9 must satisfy 
in this case is, from Eq. (|33j) . 



d9 
dt 



2C(A 3 - A) sin 9 cos 9. 



(34) 



Integrating both sides, we obtain 

9 = arctan[tan 9 exp{2C(A 3 - X)t}] . (35) 

9 = 0, 7r, and the equator (9 — tt/2) are the stationary 
points. The stability of them are determined by the sig- 
nature of A3 — A. 

The equation <fi must satisfy is, from Eqs. pip and 



^ = 2t7(A 3 - A)5'cos( 
at 



(36) 



If there is damping (a ^ 0), Eqs. (|34f and (|36|) lead to 
d6 1 



d9 a sin 9 



(37) 
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from which 



r 1 , tan (6>/2) 
6 + - In K J—t- 

ot tan(0 o /2) 



(38) 



is obtained. On the other hand, if there is no damping 
(a = 0), 9 is constant and 



(j) = 0o + 2r?(A 3 - X)St cos 9 



(39) 



is obtained. 



2. For a = 

Here we consider a case in which only the precession 
terms are present in Eqs. (|31[) - ((33]) and the three eigen- 
values arc different from each other. In this case the 
equation of motion is of the same form as that for a 
freely rotating classical rigid body, known as the Euler's 
top. The equation governing the dynamics of the top is 
called the Euler's equation [19|, whose exact solution has 
been provided by Zon and J. Schofield.poj According to 
them, we put the solution into the form, 



S x = S xm cn(uj p t + e, k), 
S y = Sy m sn(upt + e,k), 
S z = S zm dn(oj p t + £,&)) 



(40) 
(41) 
(42) 



where en, sn, and dn are the Jacobi elliptic functions. [2 1| 
The order of the eigenvalues assumed above, Eq. (|3U)) . is 
needed for this ansatz. The exact solution of the equation 
of motion is obtained by setting 



Sym = -(sgnSU) 
S zm = (sgn5 z0 ) 



1^x0 


1 A2 — A3 

Ai — A3 


o2 


\lx 2 


-^3 


o2 

y0 



Kl S*+S* , (45) 



a 3 - Ax y° 

uj p = [sgn(A 2 - A 3 )](sgnS , z0 )277 



'[(A 2 - Ar)^ + (A 3 - Ax)52 ](A 3 - A 2 ) (46) 

5; °y0 Ai — A 2 



J y0 ^ 'U^OJ 

(Ai - A 3 )5^p + (A 2 - Xs)Syo Ai - A 2 
(A 2 -A 1 )^ + (A 3 -A 1 )^ 2 ' A3-A2 



(47) 



and e = F(S y o/ S ym ; k), where F is the incomplete ellip- 
tic integral of the first kind, pi"! 



F(x;k) 



dt 



^(l-kH)(l-t 2 ) 



(48) 



The solution, Eqs. (J40J) - (|42[). indicates that, if So is 
along the x- or y- or z-axis, the spin does not change. 
3. For-q = 

Here we consider a case in which only the damping 
terms are present in Eqs. (|3"Tj) - (j3"3"]l . 

The equation <f> must satisfy is, from Eqs. ([31]) and 



dt 



2£(Ai — A 2 )cos</>sin< 



(49) 



Integrating both sides, we obtain 

4> = arctan[tan0o cxp{2C(Ai — A 2 )£}]. (50) 
Eq. (|3"3"|) thus reads 

Ai + A 2 e 4 ^ Al - A2 ) t tan 2 o " 



d6 
dt 



A, 



l + e 4C(Ai-A 2 )t tan 2 Q 

Making use of an integral formula 



Csin26». (51) 



p + qe 
1 + ae b 



dt= [ )-j-ln|l + ae ht |+pi + c, 



(52) 



where c is an integral constant, leads to 



9 = arctan 



tanfc>oe 



2C( 



Xs _ Al)t / l + e 4C(Ai-A a)ttan 2 0Q 



1 + tan^ 



(53) 



IV. CONCLUSIONS 

In conclusion, we rewrote the modified LLG equation 
applicable to both classical and quantum mechanical de- 
scription of the dynamics of spins by paying attention 
to their microscopic arrangement. The rewritten equa- 
tion of motion was demonstrated to consist of the three 
kinds of terms, each of which can be solved exactly. Their 
solutions were explicitly provided, suitable for the STD 
method. The abnormal damping term, which originates 
in the enhanced damping tensor, cannot be expressed as 
a special case of the ordinary damping term. The present 
work will help one to develop a simulation code for spin 
dynamics applicable to systems in various and complex 
situations. 
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